1 Data set up

# manon_design<- design
# manon_outcomes <- outcomes

# import data ####

# elk_rdp_ASV_iQlr.ps   <-  readRDS(file = "r_output/intermediate_rds_files/elk_rdp_ASV_iQlr.ps.rds")
# elk_rdp_genus_iQlr.ps <-  readRDS(file = "r_output/intermediate_rds_files/elk_rdp_genus_iQlr.ps.rds")
# elk_rdp_spp_iQlr.ps   <-  readRDS(file = "r_output/intermediate_rds_files/elk_rdp_spp_iQlr.ps.rds")

#CoDa conforming CLR (dont subset further)
# elk_rdp_ASV_CLR.ps    <-  readRDS(file = "r_output/intermediate_rds_files/elk_rdp_ASV_CLR.ps.rds")
# elk_rdp_day0_CLR.ps   <-  readRDS(file = "r_output/intermediate_rds_files/elk_rdp_day0_CLR.ps.rds")
# elk_rdp_genus_CLR.ps  <-  readRDS(file = "r_output/intermediate_rds_files/elk_rdp_genus_CLR.ps.rds")
elk_rdp_spp_CLR.ps    <-  readRDS(file = "r_output/intermediate_rds_files/elk_rdp_spp_CLR.ps.rds")

# elk_rdp_genus_CLR_no14.ps<- phyloseq::subset_samples(elk_rdp_genus_CLR.ps, Day != 14)


# df1.CLRday0 <- as.data.frame(elk_rdp_day0_CLR.ps@otu_table)
# df2.CLRasv <- as.data.frame(elk_rdp_ASV_CLR.ps@otu_table)
# df3.CLRgen <- as.data.frame(elk_rdp_genus_CLR.ps@otu_table)
# df3.CLRgen_no14 <- as.data.frame(elk_rdp_genus_CLR_no14.ps@otu_table)
df4.CLRspp <- as.data.frame(elk_rdp_spp_CLR.ps@otu_table)

# df1.iQlrASV <- as.data.frame(elk_rdp_ASV_iQlr.ps@otu_table)
# df2.iQlrGen  <- as.data.frame(elk_rdp_genus_iQlr.ps@otu_table)
# df3.iQlrSpp  <- as.data.frame(elk_rdp_spp_iQlr.ps@otu_table)

###
###      ID names Volunteer Sampling Tube Time
# 10111   1 10111        10        1    1    1
# 10112   2 10112        10        1    1    2
design  <- elk_rdp_spp_CLR.ps@sam_data
# design  <- elk_rdp_genus_CLR_no14.ps@sam_data
design2 <- data.frame("Day" = as.factor(design$Day),
                                            "Time" = as.numeric(design$Day),
                                          "Elk" = as.factor(design$Elk),
                                            "Rep" = as.factor(design$Rep), 
                                            "ID" = design$Description, 
                                            row.names = design$Description)

design2 <- design2 %>% mutate(names = ID) %>%
                      arrange(Day, Elk, Rep, names) %>%
                        mutate(ID = seq_along(ID)) %>%
                        mutate(ID=factor(ID, levels=ID))
rownames(design2) <- design2$names
design <- design2
n <- nrow(design)
n
## [1] 60
pander("table Elk")

table Elk

table(design$Elk)
## 
##  1  2  3  4 
## 15 15 15 15
pander("table Elk * Rep")

table Elk * Rep

table(design$Elk, design$Rep)
##    
##     1 2 3
##   1 5 5 5
##   2 5 5 5
##   3 5 5 5
##   4 5 5 5
pander("table Day * Elk")

table Day * Elk

table(design$Day, design$Elk)
##     
##      1 2 3 4
##   0  3 3 3 3
##   1  3 3 3 3
##   3  3 3 3 3
##   7  3 3 3 3
##   14 3 3 3 3

starting point can be a number of datasets

# pdf(file.path(fig_path, "HSD_outcomes.pdf"), width = 8, height = 5,
#     pointsize = 16)

outcomes_step1  <- df4.CLRspp

# outcomes_step1  <- df2.CLRasv

# outcomes_step1  <- df3.CLRgen

# outcomes_step1  <- df3.CLRgen_no14

1.1 Data dimensions and frequencies

pander("dim(design)")

dim(design)

pander(paste(dim(design),collapse = " x "))

60 x 6

pander("dim(outcomes)")

dim(outcomes)

pander(paste(dim(outcomes_step1),collapse = " x "))

60 x 162

pander("table(design$Tube, design$Sampling, design$Time)")

table(design\(Tube, design\)Sampling, design$Time)

pander(table(design$Rep, design$Elk, design$Day))
0 1 3 7 14
1 1 1 1 1 1 1
2 1 1 1 1 1
3 1 1 1 1 1
4 1 1 1 1 1
2 1 1 1 1 1 1
2 1 1 1 1 1
3 1 1 1 1 1
4 1 1 1 1 1
3 1 1 1 1 1 1
2 1 1 1 1 1
3 1 1 1 1 1
4 1 1 1 1 1
spectmat <- outcomes_step1

Elk  <- design$Elk
Rep <- design$Rep
Day <- design$Day

1.2 Hotelling T^2

model = mdatools::pca(spectmat, scale = FALSE,
            info = 'Simple PCA model', lim.type = "jm")

ncomp <- 3

Qlim <- model$Qlim
T2lim <- model$T2lim
rownames(Qlim)[2] <- rownames(T2lim)[2] <- "Out_limit"

plot(c(1:5),c(1:5))

# In case of PCA the critical limits are just shown 
# on residual plot as lines and can be used for detection 
# of extreme objects (solid line) and outliers (dashed line).

plot_hotelling <- function(){
  xlim <- range(model$calres$T2[,ncomp])
  xlim[1] <- xlim[1]*0.9
  xlim[2] <- xlim[2]*1.1
  ylim <- range(model$calres$Q[,ncomp])
  ylim[1] <- ylim[1]*0.9
  ylim[2] <- ylim[2]*2
  plot(model$calres$T2[,ncomp], model$calres$Q[,ncomp], 
     main = "Diagnostic plot for score and residual outliers", xlab = "Hotelling T2 distance", 
     ylab ="Squared residual distance", pch = 16, xlim = xlim, ylim = ylim)
  abline(h=Qlim["Out_limit",ncomp], v=T2lim["Out_limit",ncomp], lty =3)
  legend("topright", legend = "Outlier limit", lty = 3)
  index1 <- which(model$calres$T2[,ncomp]>=T2lim["Out_limit",ncomp])
  index2 <- which(model$calres$Q[,ncomp]>=Qlim["Out_limit",ncomp])
  index_ho <- unique(c(index1,index2))
  # text(x = model$calres$T2[index_ho,ncomp], y=model$calres$Q[index_ho,ncomp], 
       # labels = names(model$calres$T2[index_ho,ncomp]), pos = c(1,2,3,4))
}


# pdf(file.path(fig_path,"HSD_hotelling_CLR_spp.pdf"), height = 6, width = 6,
#     pointsize = 12)
plot_hotelling()

# dev.off()

2 Distribution of the raw outcome variables

2.1 spectmat

pander("spectmat")
par(mfrow=c(1,2))
for (i in 1:ncol(spectmat)){
  d = density(spectmat[,i])
  plot(d)
}

3 Dimension reduction by PCA

Scree plot

###################################################
# Dimension reduction by PCA
###################################################

# ===== PCA ===== #

res_pca <- SVDforPCA(x = spectmat)

### screePlot
df <- data.frame(PC =as.character(1:length(res_pca$var)),
                 var = res_pca$var)
df <- df[c(1:9),]
screePlot <- ggplot(df, aes(y=0,yend=var,x=PC,
                          xend=PC))+ geom_segment()+
          labs(title= "Scree \nplot",
                        x = "PC", y="% var")+theme_classic()
screePlot
pander("var explained per PC")

var explained per PC

pander(round(res_pca$var[1:8],2))
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
14.26 12.09 9.39 7.82 4.08 3.37 3.16 2.73

Score plots

### score plot
score1 <- DrawScores(res_pca, axes=c(1,2),drawNames = FALSE, 
           main  = "PCA scores plot", 
           pch = Day,
           color= Elk, size = 3) #+
  # coord_cartesian(xlim = c(-20, 20), ylim = c(-20, 20), expand = c(0,0))

score1.2 <- DrawScores(res_pca, axes=c(3,4),drawNames = FALSE, 
           main  = "PCA scores plot", 
           pch = Day,
           color= Elk, size = 3)#+
  # coord_cartesian(xlim = c(-20, 20), ylim = c(-20, 20), expand = c(0,0))




col <- gg_color_hue(2)
score2 <- DrawScores(res_pca, axes=c(1,2),drawNames = FALSE, 
           main  = "Scores plot", 
           color = Elk,
            size = 3)#+
  # coord_cartesian(xlim = c(-20, 20), ylim = c(-20, 20), expand = c(0,0))



# pdf(file.path(fig_path,"ScoresSpectmat.pdf"), height = 3.2, width = 9.4)

plot_grid(score1,score2, align = "none", nrow = 1)

# dev.off()

# score1
# score2

Score plot grid

score_rep1 <- DrawScores(res_pca, color = design$Rep, main="Aitchison dist. PC1 & PC2", size = 3, axes=c(1,2), drawNames = F)+
    coord_cartesian(xlim = c(-20, 20), ylim = c(-20, 20), expand = c(0,0))

score_rep2 <- DrawScores(res_pca, color = design$Rep, main="PC3 & PC4", size = 3, axes=c(3,4), drawNames = F)+
    coord_cartesian(xlim = c(-20, 20), ylim = c(-20, 20), expand = c(0,0))

score_elk1 <- DrawScores(res_pca, color = design$Elk, size = 3, axes=c(1,2), drawNames = FALSE)+
    coord_cartesian(xlim = c(-20, 20), ylim = c(-20, 20), expand = c(0,0))

score_elk2 <- DrawScores(res_pca, color = design$Elk, size = 3, axes=c(3,4), drawNames = FALSE)+
    coord_cartesian(xlim = c(-20, 20), ylim = c(-20, 20), expand = c(0,0))

score_Day1 <- DrawScores(res_pca, color = design$Day, size = 3, axes=c(1,2), drawNames = FALSE)+
    coord_cartesian(xlim = c(-20, 20), ylim = c(-20, 20), expand = c(0,0))

score_Day2 <- DrawScores(res_pca, color = design$Day, size = 3, axes=c(3,4), drawNames = FALSE)+
    coord_cartesian(xlim = c(-20, 20), ylim = c(-20, 20), expand = c(0,0))


# pdf(file.path(fig_path,"ScoresSpectmatPC1-4.pdf"), height = 12, width = 10)
plot_grid(score_rep1, score_rep2, score_elk1, score_elk2, score_Day1, score_Day2,
                    align = "none", nrow = 3)

# dev.off()

3.0.1 Graphs

## Article journal graph

loadPlot <- DrawLoadings(res_pca, type.obj = "PCA", 
                         main = "Loadings plots",
                                     axes = c(1, 2), loadingstype = "s",
                                                 ang = "90", ylab = c(""),xaxis_size = 8,
                                     xaxis_type="character",nxaxis = 20,
                                 xlab = "ASV")[[1]]

loadPlot2 <- DrawLoadings(res_pca, type.obj = "PCA", 
                         main = "Loadings plots", 
                                     axes = c(1:4), loadingstype = "s",
                                                 ang = "90", ylab = c(""),xaxis_size = 8,
                                     xaxis_type="character",nxaxis = 20,
                                 xlab = "ASV")[[1]]

PCA_plots <- ggarrange(screePlot, score1, loadPlot, 
                ncol = 3, nrow = 1,
                widths  = c(0.3, 1,  0.85))
# PCA_plots
# ggsave(PCA_plots,
#          filename = file.path(fig_path,paste0(Vers,
#            "HSD_PCA_outcomes_CLR_spp.pdf")),
#          height = 6, width = 14, scale=0.65)

score_plots_comb <- ggarrange(score1 + theme(legend.position="none"), 
                                                            score1.2+ theme(legend.position="none"), 
                ncol = 1, nrow = 2, common.legend = F,
                                                            legend = NULL)
# Extract the legend. Returns a gtable
leg <- get_legend(score1)
# Convert to a ggplot and print
leg <- as_ggplot(leg)
leg <- leg+ theme(plot.margin = unit(c(30,30,30,30), "points"))


p1 <- ggplot() + theme_void()
scree_plots_comb <- ggarrange(leg, screePlot,
                ncol = 1, nrow = 2, common.legend = F)

PCA_plots2 <- ggarrange(scree_plots_comb, score_plots_comb, loadPlot2, 
                ncol = 3, nrow = 1,
                widths  = c(0.35, 1,  0.85),
                                heights = c(0.5,1,0.95))

PCA_plots2

# ggsave(PCA_plots2,
#          filename = file.path(fig_path,paste0(Vers,
#            "HSD_PCA_outcomesPC1to4.pdf")),
#          height = 10, width = 14, scale=0.65)


score2 <- score2+
  theme(legend.key.height = unit(0.5, "cm"))

p <- ggarrange(screePlot, loadPlot,
                ncol = 2, nrow = 1,
                widths  = c(0.3,0.95))

PCA_plots3 <- ggarrange(p, score1,
                ncol = 1, nrow = 2)

# PCA_plots3

# ggsave(file.path(fig_path,paste0(Vers,"HSD_PCA_outcomes_vertical.pdf")),
#        plot = PCA_plots,
#        height = 10, width = 6.5, scale=0.75)


# score1 <- score1+
#   theme(legend.key.height = unit(0.5, "cm"))
# score1.2 <- score1.2+
#   theme(legend.key.height = unit(0.5, "cm"))
# 
# p <- ggarrange(screePlot, loadPlot2,
#                 ncol = 2, nrow = 1,
#                 widths  = c(0.3,0.95))
# 
# PCA_plots <- ggarrange(p, score1, score1.2,
#                 ncol = 1, nrow = 3, heights = c(1,1.3,1.3),
#                                               common.legend = T, legend = "bottom")
# 
# PCA_plots

3.1 Apply first step of LiMMPCA (PCA)

3.1.1 transform outcomes??

# transform outcomes??
transfo = FALSE
if (transfo){
  outcomes <- log(outcomes + 0.01)
}
# apply first step of LiMMPCA with dimreducPCA()
res_pca2 <- dimreducPCA(data = as.matrix(outcomes_step1), pcvar = 95)


spectra_PCA_scores <- res_pca2$pca_scores
spectra_PCA_loadings <- res_pca2$pca_loadings


m <- nPC <- dim(spectra_PCA_scores)[2]
pander("Number of PCs kept:")

Number of PCs kept:

m
## [1] 40
# balanced ??
balanced=FALSE

if (balanced){
  spectra_PCA_scores_bal <-  spectra_PCA_scores[rownames(spectra_PCA_scores) %in% 
                       rownames(design_balanced), ]
  design <- design_balanced
  spectra_PCA_scores <- spectra_PCA_scores_bal
}

4 Parallel mixed modelling

###################################################
# Parallel mixed modelling
###################################################

### define the input arguments for the parallel mixed modelling

# formula
# form <- "~ Time + Tube + (1|Volunteer) + (1|Volunteer:Sampling) "
# form <- "~ Day + (1|Rep) + (1|Elk)" #ideal
form <- "~ Day * Elk + (1|Rep)"
# form <- "~ Day * Elk + (1|Rep) + (1|Elk) + (1|Rep:Elk)"
# form <- "~ Day +  (1|Elk) + (1|Elk:Day) + (1|Rep) + (1|Elk:Rep) + (1|Day:Rep)" #ideal
# form <- "~ Day + (1|Elk/Rep)" # NO NEED for this, it explains the same variation as above
# form <- "~ Time + (1|Elk) + (1|Rep) " #cont variables do not work
# multivariate response matrix
outcomes_step2 <- spectra_PCA_scores

# REML
REML <- TRUE

### run parlmer: lmer in parrallel on each response vector
res.parlmer <- parlmer(design, outcomes_step2, form, REML)
MM_full <- res.parlmer


# save the results
RanModMatlist <- res.parlmer$RanModMatlist
FixedModMatlist <- res.parlmer$FixedModMatlist

# Residuals sd error
Res_std_error_PC <- sapply(MM_full$merMod_obj, sigma)

# Residuals
Residuals_PC <- sapply(MM_full$merMod_obj, residuals)


### ranef_PC: Extract the modes of the random effects
ranef_PC <- sapply(MM_full$merMod_obj, function(x) unlist(ranef(x)))


# recover accurate rownames and colnames of ranef_PC
list_rownam <- lapply(ranef(MM_full$merMod_obj[[1]]), rownames)

colnam <- paste0(names(ranef(MM_full$merMod_obj[[1]])), lapply(ranef(MM_full$merMod_obj[[1]]), rownames))

ranef_PC <- sapply(MM_full$merMod_obj, function(x) unlist(ranef(x)))
names(list_rownam) <- gsub("[^A-z]", "", names(list_rownam))

colnam <- c()
for (i in 1:length(list_rownam)){
  colnam <- c(colnam, paste0(names(list_rownam)[i], list_rownam[[i]]))
}
rownames(ranef_PC) <- colnam
rownames(ranef_PC) <- gsub("..(Intercept))", "", rownames(ranef_PC))


### fixef: Extract fixed-effects estimates
fixef_PC <- sapply(MM_full$merMod_obj, fixef)


### all fixed estimates and random predictions
cof_PC <- rbind(fixef_PC, ranef_PC)

### Extract Variance and Correlation Components
varcor_random_full <- sapply(MM_full$merMod_obj, VarCorr) # var

dat <- as.numeric(rbind(varcor_random_full))
varcor_random_full_mat <- matrix(dat, nrow=length(varcor_random_full)[1],
                                 dimnames = list(names(varcor_random_full)))


fixNames <- MM_full$fixNames # names of fixed effects

ranNames <- MM_full$ranNames # names of random effects

# estimated Std.Dev. of fixed effects parameters
varcor_fixed_full <- sapply(MM_full$merMod_obj, 
                            function(x) sqrt(diag(vcov(x)))) # Std.Dev.


rownames(varcor_fixed_full) <- rownames(vcov(MM_full$merMod_obj[[1]]))

5 Normality of the model residuals

# pdf(file.path(out_path, "HSD_residuals.pdf"), height = 15, width = 10)
par(mfrow=c(8,4), mar=c(2,2,2,2))
ncol(Residuals_PC) #51
## [1] 40
for (i in 1:15){
  qqPlot(Residuals_PC[,i], main = paste0("qqplot - PC",i), 
         ylab="Residuals quantiles")
  d = density(Residuals_PC[,i])
  hist(Residuals_PC[,i], main = paste0("Histogram - PC",i), xlab="",
       breaks = 20)
}
# dev.off()

6 Effect matrix computation

###################################################
# Effect matrix computation
###################################################

# fixed effects + intercept -----------------------------
dim1FixedModMad <- sapply(FixedModMatlist, function(x) dim(x)[2])
names_FixedEffects <- names(FixedModMatlist)

shortFixedNames <- names_FixedEffects
Xmat <- do.call(cbind, FixedModMatlist)

# colnames(Xmat) %in% rownames(fixef_PC)
Xmat <- Xmat[,rownames(fixef_PC)] # reorder colnames of Xmat

index <- cumsum(dim1FixedModMad)
k <- 1
Mfix <- vector("list", length=length(shortFixedNames))
Mfix_PC <- vector("list", length=length(shortFixedNames))
names(Mfix) <- names(Mfix_PC) <-  shortFixedNames
for (i in 1:length(shortFixedNames)){
  XMfix = Xmat
  XMfix[,-c(k:index[i])] = 0
  Mfix_PC[[i]] = XMfix %*% fixef_PC
  Mfix[[i]] <- Mfix_PC[[i]]%*%t(spectra_PCA_loadings)
  k <-  index[i] + 1
  dimnames(Mfix[[i]]) <- dimnames(spectmat)
  rownames(Mfix_PC[[i]]) <- rownames(spectmat)
}

M0 <- Mfix$Intercept

Mfix <- Mfix[-which(names(Mfix)=="(Intercept)")]
Mfix_PC <- Mfix_PC[-which(names(Mfix_PC)=="(Intercept)")]


# random effects -----------------------------
dim1RandModMad <- sapply(RanModMatlist, function(x) dim(x)[2])
names_randomEffects <- names(RanModMatlist)

shortRandNames <- gsub("[^A-z]", "", names_randomEffects)

Zmat <- do.call(cbind, RanModMatlist)
colnames(Zmat) <- paste0(rep(shortRandNames,
                             dim1RandModMad),
                         colnames(Zmat))

# colnames(Zmat) %in% rownames(ranef_PC)

Zmat <- Zmat[,rownames(ranef_PC)] # reorder colnames of Zmat

index <- cumsum(dim1RandModMad)
k <- 1
Mrand_PC <- vector("list", length=length(ranNames))
Mrand <- vector("list", length=length(ranNames))
names(Mrand) <- names(Mrand_PC) <- ranNames
for (i in 1:length(shortRandNames)){
  XMrand = Zmat
  XMrand[,-c(k:index[i])] = 0
  Mrand_PC[[i]] = XMrand %*% ranef_PC
  Mrand[[i]] <- Mrand_PC[[i]]%*%t(spectra_PCA_loadings)
  k <-  index[i] + 1
  dimnames(Mrand[[i]]) <- dimnames(spectmat)
  rownames(Mrand_PC[[i]]) <- rownames(spectmat)
}

Mres_PC <- Residuals_PC
Mres <- Mres_PC%*%t(spectra_PCA_loadings)

7 PCA on the effect matrices

7.1 PCA on the pure effect matrices

7.1.1 Scree plots and pure loadings

Note that only the loadings are backtransformed

# ======================================
# Scree plots and pure loadings ========
# ======================================

# Time
#####################
res_pca <- SVDforPCA(Mfix_PC$Day)

df <- data.frame(PC = as.character(1:7),var = res_pca$var[1:7])

screeplotTime <- ggplot(df, aes(y=0,yend=var,x=PC,
                          xend=PC))+ geom_segment()+
                   labs(title= "Scree plot",
                        x = "PC", y="% var")+theme_classic()
    
pander("screeplot Time")       

screeplot Time

screeplotTime

backtransf_load <- t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))

colnames(backtransf_load) <- paste0("L", 
                          1:ncol(backtransf_load), 
                          " (",round(res_pca$var,2),"%)")

loadTime <- LinePlot(t(backtransf_load), rows=c(1:4), 
                                         nxaxis = 20,ang = "90", xaxis_type = "character",
                     main = "PCA loadings",  
                     xlab = "asv", type = "s")[[1]]

pander("loadings Time") 

loadings Time

loadTime

# Elk
#####################
res_pca <- SVDforPCA(Mfix_PC$Elk)

df <- data.frame(PC = as.character(1:7),var = res_pca$var[1:7])

screeplotElk <- ggplot(df, aes(y=0,yend=var,x=PC,
                          xend=PC))+ geom_segment()+
                   labs(title= "Scree plot",
                        x = "PC", y="% var")+theme_classic()
    
pander("screeplot Elk")       

screeplot Elk

screeplotElk

backtransf_load <- t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))

colnames(backtransf_load) <- paste0("L", 
                          1:ncol(backtransf_load), 
                          " (",round(res_pca$var,2),"%)")

loadElk <- LinePlot(t(backtransf_load), rows=c(1:3), 
                                         nxaxis = 20,ang = "90", xaxis_type = "character",
                     main = "PCA loadings",  
                     xlab = "asv", type = "s")[[1]]

pander("loadings Elk") 

loadings Elk

loadElk

# Sampling
#####################
res_pca <- SVDforPCA(Mrand_PC$Rep)

df <- data.frame(PC = as.character(1:7),var = res_pca$var[1:7])
screeplotSampling <- ggplot(df, aes(y=0,yend=var,x=PC,
                          xend=PC))+ geom_segment()+
                   labs(title= "Scree plot",
                        x = "PC", y="% var")+theme_classic()
   
pander("screeplot Sampling")         

screeplot Sampling

screeplotSampling

backtransf_load <- t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))
colnames(backtransf_load) <- paste0("L", 
                          1:ncol(backtransf_load), 
                          " (",round(res_pca$var,2),"%)")


loadSampling <- LinePlot(t(backtransf_load), rows=c(1,2), 
                                         nxaxis = 20,ang = "90", xaxis_type = "character",
                     main = "PCA loadings",  
                     xlab = "asv", type = "s")[[1]]

pander("loadings Sampling")

loadings Sampling

loadSampling

# # Volunteer
# #####################
# res_pca <- SVDforPCA(Mrand_PC$Elk)
# 
# df <- data.frame(PC = as.character(1:7),var = res_pca$var[1:7])
# screeplotVolunteer <- ggplot(df, aes(y=0,yend=var,x=PC,
#                           xend=PC))+ geom_segment()+
#                    labs(title= "Scree plot",
#                         x = "PC", y="% var")+theme_classic()
# pander("screeplot Volunteer")  
# screeplotVolunteer
# 
# backtransf_load <- t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))
# colnames(backtransf_load) <- paste0("L", 
#                           1:ncol(backtransf_load), 
#                           " (",round(res_pca$var,2),"%)")
# 
# backtransf_load[,2] <- backtransf_load[,2] *-1
# 
# loadVolunteer <- LinePlot(t(backtransf_load), rows=c(1,2,3), 
#                                       nxaxis = 20,ang = "90", xaxis_type = "character",
#                      main = "PCA loadings",  
#                      xlab = "asv", type = "s")[[1]]
# pander("loadings Volunteer")
# loadVolunteer

7.2 PCA on the augmented effect matrices

7.2.1 addSegments function

addSegments <- function(group, data, pch=16, main = NULL, 
                        col = rainbow(n = length(unique(group))),
                        ...) {
  
  group <- as.factor(group)
  
  plot(data$x, data$y, col=col[group], pch=pch, main=main, ...)
  
  group <- as.factor(group)
  xcent <- tapply(data[,1], group, FUN=mean)
  ycent <- tapply(data[,2], group, FUN=mean)
  centers <- data.frame(xcent=xcent, ycent=ycent)

  mapply(FUN = points, x = centers$xcent, 
         y = centers$ycent, col = col, 
         MoreArgs = list(pch=20, cex=0.7))

  submatrices <- split(x=data, f=group)
  for (i in 1:nlevels(group)){
    mapply(FUN = segments, x1 = submatrices[[i]]$x, 
           y1 = submatrices[[i]]$y, col = col[i],
         MoreArgs = list(x0 = centers$xcent[i], 
                         y0 = centers$ycent[i]))
  }
  
}

7.2.2 Augmented (with ANOVA degrees of freedom) scores plots (ASCA-E)

# ========================
# augmented scores =======
# ========================
# * Tube: Mu + E
# * Time: Mt + E
# * Volunteer: Mv + Ms
# * Sampling: Ms + E

# correction factors ++++++
# a <- nlevels(design$Tube)
# b <- nlevels(design$Time)
# c <- nlevels(design$Volunteer)
# d <- nlevels(design$Sampling)

b <- nlevels(design$Day)
c <- nlevels(design$Elk)
d <- nlevels(design$Rep)

#################################
###  Augmented time effect
#################################

pander("Augmented Time effect")

Augmented Time effect

res_pca <- SVDforPCA(Mfix_PC$Day)
df1 <- (b-1)
df2 <- (b*c*d-b-c*d+2)
Fstat <- qf(.95, df1=df1, df2=df2) 
coef <- sqrt(Fstat*df1/df2)

res_pca$scores[,1:2] <- (Mfix_PC$Day + Mres_PC*coef)%*%res_pca$loadings[,1:2]

col <- c("1"=violetred, "2" = darkblue)
pch <- c("1"= 2, "2"= 4)
# index <- c(84, 85, 87, 88)
# PC1lab <- paste0("PC", 1, " (", round(res_pca$var[1], 2),"%)")
df <- as.data.frame(res_pca$scores)

Time_scores <- ggplot(aes(y=PC1, x= Day, color=Day, 
                          shape = Elk),data=df)+
  geom_jitter() +
  ggtitle(label = "PCA Scores")+
  xlab("Day")+  #ylim(low=-2.1, high=2.1) +
  geom_hline(yintercept=0, linetype=2, color="gray", size=0.7)+
  theme(plot.title = element_text(face="plain"),axis.text.x = element_text(size=10)) +
    theme_bw()
   # annotate("text", y = res_pca$scores[index,1]+ 0.1*c(1,1,1,1), 
   #       x = index, label = rownames(res_pca$scores)[index])
Time_scores

#################################
###  Augmented Elk effect
#################################

pander("Augmented Elk effect")

Augmented Elk effect

res_pca <- SVDforPCA(Mfix_PC$Elk)
df1 <- (c-1)
df2 <- (b*c*d-b-c*d+2)
Fstat <- qf(.95, df1=df1, df2=df2) 
coef <- sqrt(Fstat*df1/df2)

res_pca$scores[,1:2] <- (Mfix_PC$Elk + Mres_PC*coef)%*%res_pca$loadings[,1:2]

col <- c("1"=violetred, "2" = darkblue)
pch <- c("1"= 2, "2"= 4)
# index <- c(84, 85, 87, 88)

# PC1lab <- paste0("PC", 1, " (", round(res_pca$var[1], 2),"%)")

df <- as.data.frame(res_pca$scores)

Elk_scores <- ggplot(aes(y=PC1, x= Day, color=Elk, group = Elk,
                          shape = Day),data=df)+
  geom_jitter() +
    # geom_point()+
    labs(title = "PCA Scores", 
                x = "Day", y = "PC1") +
  geom_hline(yintercept=0, linetype=2, color="gray", size=0.7)+
  theme(plot.title = element_text(face="plain"),axis.text.x = element_text(size=10))+
    theme_bw()
   # annotate("text", y = res_pca$scores[index,1]+ 0.1*c(1,1,1,1), 
   #       x = index, label = rownames(res_pca$scores)[index])
Elk_scores

#plot 2
xlab <- paste0("PC", 1, " (", round(res_pca$var[1], 2),"%)")
ylab <- paste0("PC", 2, " (", round(res_pca$var[2], 2), "%)")

pch <- c("1"= 1, "2"= 17, "3"=8)

col <- c(darkblue, "dodgerblue", turquoise, violetred, 
         limegreen,"orangered", "forestgreen", "gold", 
         "goldenrod4", "orange1", "maroon1",
         "red3")

VolSamp <- as.factor(paste0("elk_",Elk,".day_", Day))
Vol <- as.factor(substr(VolSamp,1,5))
# pch <- c("1"= 2, "2"= 8)
names(col) <- levels(Elk)

Elk_scores2 <- ggplot(df, aes(PC1, PC2, group = VolSamp, 
               shape = Day)) +   
  geom_point(size = 2, aes(colour = Elk)) +
  geom_line(aes(colour = Elk)) +
  # scale_shape_manual(name = "Day", values=c(8, 0,17, 1, 2),
  #                          guide=guide_legend(order=1, shape = 1)) + 
  # coord_cartesian(xlim = c(-5, 5), ylim = c(-10, 10))   + 
    
  ggplot2::labs(title = "PCA Scores", 
                x = xlab, y = ylab) +
  ggplot2::geom_vline(xintercept = 0, size = 0.1) + 
  ggplot2::geom_hline(yintercept = 0, size = 0.1) + 
  ggplot2::theme_bw() +
      ggplot2::theme(panel.grid.major = 
                       ggplot2::element_line(color = "gray60",
                       size = 0.2), 
                     panel.grid.minor = ggplot2::element_blank(),
                     panel.background = 
                       ggplot2::element_rect(fill = "gray98")) +   theme(legend.text=element_text(size=7),
                    legend.key.height=unit(0.55,"line"))
# Elk_scores2
# #################################
# ### Augmented Volunteer effect
# #################################
# 
# pander("Augmented Volunteer effect")
# 
# res_pca <- SVDforPCA(Mrand_PC$Elk)
# df1 <- (c-1)
# df2 <- (c*d-c)
# Fstat <- qf(.95, df1=df1, df2=df2) 
# coef <- sqrt(Fstat*df1/df2)
# 
# res_pca$scores[,1:2] <- (Mrand_PC$Elk+
#                            Mres_PC*coef)%*%
#   res_pca$loadings[,1:2]
# res_pca$scores[,2] <- res_pca$scores[,2]*-1
# 
# res_pca$scores <- round(res_pca$scores, 6)
# Volunteer_scores <- DrawScores(res_pca, 
#                main = "PCA Scores", 
#            color = Elk, pch=Day, 
#            drawNames = FALSE, drawPolygon  = TRUE, 
#            noLegend = FALSE, size=2)+
#         coord_cartesian(xlim = c(-25, 25), 
#                         ylim = c(-25, 25))+  
#   theme(legend.text=element_text(size=10),
#         plot.title = element_text(face="plain"),
#                     legend.key.height=unit(0.55,"line"))
# 
# Volunteer_scores
# #################################
# ### Augmented Day effect
# #################################
# 
# pander("Augmented Volunteer effect")
# 
# res_pca <- SVDforPCA(Mrand_PC$Day)
# df1 <- (b-1)
# df2 <- (b*c*d-b-c*d+2)
# Fstat <- qf(.95, df1=df1, df2=df2) 
# coef <- sqrt(Fstat*df1/df2)
# 
# res_pca$scores[,1:2] <- (Mrand_PC$Day+
#                            Mres_PC*coef)%*%
#   res_pca$loadings[,1:2]
# res_pca$scores[,2] <- res_pca$scores[,2]*-1
# 
# res_pca$scores <- round(res_pca$scores, 6)
# Volunteer_scores <- DrawScores(res_pca, 
#                main = "PCA Scores", 
#            color = Elk, pch=Day, 
#            drawNames = FALSE, drawPolygon  = TRUE, 
#            noLegend = FALSE, size=2)+
#         coord_cartesian(xlim = c(-25, 25), 
#                         ylim = c(-25, 25))+  
#   theme(legend.text=element_text(size=10),
#         plot.title = element_text(face="plain"),
#                     legend.key.height=unit(0.55,"line"))
# 
# Volunteer_scores
#################################
### Augmented Sampling effect
#################################
pander("Augmented Sampling effect")

Augmented Sampling effect

res_pca <- SVDforPCA(Mrand_PC$Rep)
df1 <- (d-1)
df2 <- (b*c*d-b-c*d+2)
Fstat <- qf(.95, df1=df1, df2=df2) 
coef <- sqrt(Fstat*df1/df2)

res_pca$scores[,1:2] <- (Mrand_PC$Rep + Mres_PC*coef)%*%res_pca$loadings[,1:2]

res_pca$scores <- round(res_pca$scores, 6)
Sampling_scores <- DrawScores(res_pca, 
               main = "PCA Scores", 
           color = Rep, pch=Elk, 
           drawNames = FALSE, drawPolygon  = TRUE, 
           noLegend = FALSE, size=2)+
        coord_cartesian(xlim = c(-25, 25), 
                        ylim = c(-25, 25))+  
  theme(legend.text=element_text(size=10),
        plot.title = element_text(face="plain"),
                    legend.key.height=unit(0.55,"line"))

pch <- c("1"= 1, "2"= 17, "3"=8)

col <- c(darkblue, "dodgerblue", turquoise, violetred, 
         limegreen,"orangered", "forestgreen", "gold", 
         "goldenrod4", "orange1", "maroon1",
         "red3")

VolSamp <- as.factor(paste0("elk_",Elk,".rep_", Rep))
Vol <- as.factor(substr(VolSamp,1,5))
pch <- c("1"= 2, "2"= 8)
names(col) <- levels(Elk)

xlab <- paste0("PC", 1, " (", round(res_pca$var[1], 2),"%)")
ylab <- paste0("PC", 2, " (", round(res_pca$var[2], 2), "%)")
    
df <- as.data.frame(res_pca$scores[,1:2])
Sampling_scores <- ggplot(df, aes(PC1, PC2, group = VolSamp, 
               shape = Rep)) +   
  geom_point(size = 2, aes(colour = Elk)) +
  geom_line(aes(colour = Elk)) +
  scale_shape_manual(name = "Sampling", values=c(8, 0,17),
                           guide=guide_legend(order=1, shape = 1)) + 
  coord_cartesian(xlim = c(-5, 5), ylim = c(-10, 10))   + 
    
  ggplot2::labs(title = "PCA Scores", 
                x = xlab, y = ylab) +
  ggplot2::geom_vline(xintercept = 0, size = 0.1) + 
  ggplot2::geom_hline(yintercept = 0, size = 0.1) + 
  ggplot2::theme_bw() +
      ggplot2::theme(panel.grid.major = 
                       ggplot2::element_line(color = "gray60",
                       size = 0.2), 
                     panel.grid.minor = ggplot2::element_blank(),
                     panel.background = 
                       ggplot2::element_rect(fill = "gray98")) +   theme(legend.text=element_text(size=7),
                    legend.key.height=unit(0.55,"line")) + theme_bw()


Sampling_scores

7.2.3 Augmented (with Effective Dimensions) scores plots (ASCA-E)

# Effective Dimensions
# source("ED/metabo_Mixed.R")
EDu <- rep((nlevels(design$Rep)-1),ncol(outcomes_step2))
EDt <- rep((nlevels(design$Day)-1),ncol(outcomes_step2))
# EDs <- ED_metabo["ED_Samp",]
# EDv <- ED_metabo["ED_Vol",]
# EDe <- n - colSums(ED_metabo)
EDe <- 20

# ========================
# augmented scores =======
# ========================
# * Tube: Mu + E
# * Time: Mt + E
# * Volunteer: Mv + Ms
# * Sampling: Ms + E


#################################
###  Augmented time effect
#################################

# pander("Augmented Time effect")
# res_pca <- SVDforPCA(Mfix_PC$Day)
# df1 <- EDt
# df2 <- EDe
# Fstat <- qf(.95, df1=df1, df2=df2) 
# coef <- sqrt(Fstat*df1/df2)
# 
# mat <- matrix(NA, ncol=ncol(outcomes_step2), nrow=nrow(outcomes_step2))
# for (i in 1:ncol(outcomes_step2)){
#   mat[,i] <- Mres_PC[,i]*coef[i]
# }
# 
# res_pca$scores[,1:2] <- (Mfix_PC$Day + mat)%*%res_pca$loadings[,1:2]
# 
# col <- c("1"=violetred, "2" = darkblue)
# pch <- c("1"= 2, "2"= 4)
# # index <- c(84, 85, 87, 88)
# 
# df <- as.data.frame(res_pca$scores)
# 
# Time_scores <- ggplot(aes(y=PC1, x= 1:nrow(df), color=Day, 
#                           shape = Rep),data=df)+
#   geom_point() +
#   ggtitle(label = "PCA Scores")+
#   xlab("Observation number")+  #ylim(low=-2.1, high=2.1) +
#   geom_hline(yintercept=0, linetype=2, color="gray", size=0.7)+
#   theme(plot.title = element_text(face="plain"),axis.text.x = element_text(size=10)) #+
#    # annotate("text", y = res_pca$scores[index,1]+ 0.1*c(1,1,1,1), 
#    #       x = index, label = rownames(res_pca$scores)[index])
# Time_scores



#################################
###  Augmented tube effect
#################################

# pander("Augmented Tube effect")
# res_pca <- SVDforPCA(Mfix_PC$Elk)
# df1 <- EDu
# df2 <- EDe
# Fstat <- qf(.95, df1=df1, df2=df2) 
# coef <- sqrt(Fstat*df1/df2)
# 
# mat <- matrix(NA, ncol=ncol(outcomes_step2), nrow=nrow(outcomes_step2))
# for (i in 1:ncol(outcomes_step2)){
#   mat[,i] <- Mres_PC[,i]*coef[i]
# }
# 
# res_pca$scores[,1:2] <- (Mfix_PC$Elk + mat)%*%res_pca$loadings[,1:2]
# 
# col <- c("1"=violetred, "2" = darkblue)
# pch <- c("1"= 2, "2"= 4)
# 
# # index <- c(80,79, 84,81,78,22,19,18,21,87)
# df <- as.data.frame(res_pca$scores)
# 
# Tube_scores <- ggplot(aes(y=PC1, x= 1:nrow(df), 
#                           color=Elk, shape = Day),
#                       data=df) +
#   geom_point() +
#   ggtitle(label = "PCA Scores")+
#   xlab("Observation number")+
#   geom_hline(yintercept=0, linetype=2, color="gray", size=0.7)+
#   theme(plot.title = element_text(face="plain"),
#         axis.text.x = element_text(size=10)) #+
#    # annotate("text", 
#    #          y = res_pca$scores[index,1]+ 0.2*
#    #            c(1,1,-1,1,-1,1, -1,-1,1,-1), 
#    #       x = index, label = rownames(res_pca$scores)[index])
# 
# Tube_scores


#################################
### Augmented Volunteer effect
#################################

# pander("Augmented Volunteer effect")
# 
# res_pca <- SVDforPCA(Mrand_PC$Volunteer)
# df1 <- EDv
# df2 <- EDs
# Fstat <- qf(.95, df1=df1, df2=df2) 
# coef <- sqrt(Fstat*df1/df2)
# 
# mat <- matrix(NA, ncol=ncol(outcomes), nrow=nrow(outcomes))
# for (i in 1:ncol(Mrand_PC$Sampling)){
#   mat[,i] <- Mrand_PC$Sampling[,i]*coef[i]
# }
# 
# res_pca$scores[,1:2] <- (Mrand_PC$Volunteer+mat)%*%
#   res_pca$loadings[,1:2]
# res_pca$scores[,2] <- res_pca$scores[,2]*-1
# 
# res_pca$scores <- round(res_pca$scores, 6)
# Volunteer_scores <- DrawScores(res_pca, 
#                main = "PCA Scores", 
#            color = Volunteer, pch=Volunteer, 
#            drawNames = FALSE, drawPolygon  = TRUE, 
#            noLegend = FALSE, size=2)+
#         coord_cartesian(xlim = c(-50, 50), 
#                         ylim = c(-17, 17))+  
#   theme(legend.text=element_text(size=10),
#         plot.title = element_text(face="plain"),
#                     legend.key.height=unit(0.55,"line"))
# 
# Volunteer_scores

#################################
### Augmented Sampling effect
#################################
# pander("Augmented Sampling effect")
# 
# res_pca <- SVDforPCA(Mrand_PC$Sampling)
# df1 <- EDs
# df2 <- EDe
# Fstat <- qf(.95, df1=df1, df2=df2) 
# coef <- sqrt(Fstat*df1/df2)
# 
# mat <- matrix(NA, ncol=ncol(outcomes), nrow=nrow(outcomes))
# for (i in 1:ncol(outcomes)){
#   mat[,i] <- Mres_PC[,i]*coef[i]
# }
# 
# res_pca$scores[,1:2] <- (Mrand_PC$Sampling + mat)%*%res_pca$loadings[,1:2]
# 
# res_pca$scores <- round(res_pca$scores, 6)
# 
# pch <- c("1"= 1, "2"= 17, "3"=8)
# 
# col <- c(darkblue, "dodgerblue", turquoise, violetred, 
#          limegreen,"orangered", "forestgreen", "gold", 
#          "goldenrod4", "orange1", "maroon1",
#          "red3")
# 
# VolSamp <- as.factor(paste0(Volunteer, Sampling))
# Vol <- as.factor(substr(VolSamp,1,2))
# pch <- c("1"= 2, "2"= 8)
# names(col) <- levels(Volunteer)
# 
# xlab <- paste0("PC", 1, " (", round(res_pca$var[1], 2),"%)")
# ylab <- paste0("PC", 2, " (", round(res_pca$var[2], 2), "%)")
#     
# df <- as.data.frame(res_pca$scores[,1:2])
# Sampling_scores <- ggplot(df, aes(PC1, PC2, group = VolSamp, 
#                shape = Sampling)) +   
#   geom_point(size = 2, aes(colour = Volunteer)) +
#   geom_line(aes(colour = Volunteer)) +
#   scale_shape_manual(name = "Sampling", values=c(8, 0,17),
#                            guide=guide_legend(order=1, shape = 1)) + 
#   coord_cartesian(xlim = c(-26, 26), ylim = c(-12, 12))   + 
#   ggplot2::labs(title = "PCA Scores", 
#                 x = xlab, y = ylab) +
#   ggplot2::geom_vline(xintercept = 0, size = 0.1) + 
#   ggplot2::geom_hline(yintercept = 0, size = 0.1) + 
#   ggplot2::theme_bw() +
#       ggplot2::theme(panel.grid.major = 
#                        ggplot2::element_line(color = "gray60",
#                        size = 0.2), 
#                      panel.grid.minor = ggplot2::element_blank(),
#                      panel.background = 
#                        ggplot2::element_rect(fill = "gray98")) +   theme(legend.text=element_text(size=10),
#                     legend.key.height=unit(0.55,"line"))
# 
# 
# Sampling_scores

Plots

### Thesis chapter output
Time_scores <- Time_scores + theme(text=element_text(size=12))
loadTime <- loadTime + theme(text=element_text(size=12))
a <- grid.arrange(screeplotTime, Time_scores,
                 loadTime, nrow=1,widths=c(0.3, 1,  1),
             top=textGrob("Time effect matrix",
                          gp=gpar(fontsize=20,font=2)))

Elk_scores <- Elk_scores + theme(text=element_text(size=12))
loadElk <- loadElk + theme(text=element_text(size=12))
b <-  grid.arrange(screeplotElk,
                  Elk_scores,loadElk, 
                  nrow=1,widths=c(0.3, 1,  1),
             top=textGrob("Elk effect matrix",
                          gp=gpar(fontsize=20,font=2)))

Sampling_scores <- Sampling_scores + theme(text=element_text(size=12))
loadSampling <- loadSampling + theme(text=element_text(size=12))
c <- grid.arrange(screeplotSampling,
               Sampling_scores,loadSampling,
              nrow=1,widths=c(0.3, 1,  1),
             top=textGrob("Sampling effect matrix",
                          gp=gpar(fontsize=20,font=2)))

## Thesis chapter output
# pdf(file.path(fig_path,paste0(Vers,
           # "HSD_Scores_Loadings_EffectMat_backtrans.pdf")),
         # height = 20, width = 14)
ggarrange(a,b,c, nrow=3, labels = c("A", "B", "C"))

# dev.off()


# plots <- ggarrange(a,b, nrow=2, labels = c("A", "B"))
# plots
# ggsave(file.path(fig_path,paste0(Vers,
           # "HSD_Scores_Loadings_EffectMat_backtrans1.pdf")), plots,
        # height = 10, width = 15, scale=0.68)

# plots <- ggarrange(c,d, nrow=2, labels = c("C", "D"))
# plots
# ggsave(file.path(fig_path,paste0(Vers,
           # "HSD_Scores_Loadings_EffectMat_backtrans2.pdf")), plots,
        # height = 10, width = 15, scale=0.68)

7.2.4 Residuals

res_pca <- SVDforPCA(Mres_PC)

df <- data.frame(PC = as.character(1:length(res_pca$var)),
                 var = res_pca$var)

df <- df[1:7,]

screeplot_resid <- ggplot(df, aes(y=0,yend=var,x=PC,
                          xend=PC))+ geom_segment() +
          labs(title= "Scree plot",
                        x = "PC", y="% var") + theme_classic()

backtransf_load <- t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))

colnames(backtransf_load) <- paste0("L", 
                          1:ncol(backtransf_load), 
                          " (",round(res_pca$var,2),"%)")

loadings_resid <- LinePlot(t(backtransf_load), rows=c(1,2), 
                     main = "PCA loadings", 
                                        xaxis_type = "character",ang = "90", nxaxis = 20,
                     xlab = "asv", type = "s")[[1]]
loadings_resid <- loadings_resid + theme(text=element_text(size=12))

# index <- c(85,  84, 80, 81, 79, 78)


res_pca <- SVDforPCA(Mres_PC)
scores_resid <- DrawScores(res_pca, drawNames = FALSE, 
                        # xaxis_type = "character",ang = "90", nxaxis = 20,
            color = Day, 
           pch = Elk, main ="PCA scores", 
           size = 2) + 
  # annotate("text", y = (res_pca$scores[index,2]+0.3 +
  #                               0.4*c( 1 , 1, 1, 1,-1,1)),
  #                        x = res_pca$scores[index,1]+0.3,
  #                        label = rownames(res_pca$scores[index,1:2]))+ 
   theme(legend.text=element_text(size=10),
         text=element_text(size=12),
         legend.key.height=unit(0.7,"line")) 


# scores_resid <- scores_resid + coord_cartesian(xlim = c(-9, 9), 
                                               # ylim = c(-9, 9))
# scores_resid


plots <- grid.arrange(screeplot_resid,
              scores_resid, loadings_resid,
              nrow=1,widths  = c(0.3, 1,  1),
             top=textGrob("Residual effect matrix",
                          gp=gpar(fontsize=20,font=2)))

# plots

# ggsave(filename = paste0(fig_path,"/HSD_Scores_Loadings_Residuals.pdf"),
             # plot = plots, device = "pdf", height = 6, width = 12, scale = 0.9)

plots combined

### Journal article output
e <- grid.arrange(screeplot_resid, 
              scores_resid, loadings_resid,
              nrow=1,widths  = c(0.3, 1,  1),
             top=textGrob("Residual effect matrix",
                          gp=gpar(fontsize=20,font=2)))

# journal article output

# pdf(file.path(fig_path,paste0(Vers,
#            "HSD_Scores_Loadings_EffectMat.pdf")),
#          height = 22, width = 14)
p <- ggarrange(a,b,c,e , nrow=4, labels = c("A", "B", "C", "D"))
p

# dev.off()

# ggsave(file.path(fig_path,paste0(Vers,
           # "HSD_Scores_Loadings_EffectMat.pdf")), plot = p,
       # height = 22, width = 14, scale = 0.8)

8 Percentage of explained variance

8.1 Method from Nakagawa and Schielzeth (2012)

####################################################
# Pourcentage de variance expliquée
####################################################

# Random effects -----------------------------------
sigma2_res = Res_std_error_PC^2 # Residual
varcor_random_full <- as.data.frame(varcor_random_full)
Var_Mrand <- rbind(varcor_random_full, 
                   sigma2_res=sigma2_res) 
Var_Mrand <- data.matrix(Var_Mrand)

ranNames <- rownames(varcor_random_full)

# fixed effects -----------------------------------
# Mfix_PC[1]
Var_Mfix <- c()
for (i in 1:length(fixNames)){
  #  variance of parameter values (population)
  Var_Mfix <- rbind(Var_Mfix, 
                    (apply(Mfix_PC[[i]], 2, var) *(n - 1) / n))
}
rownames(Var_Mfix) <- fixNames

rownames(Var_Mrand) <- c(ranNames, "Residuals")

Var_Mrand <- rbind(
                                     # Elk=Var_Mrand["Elk",],
                                     # Day=Var_Mrand["Day",],
                                     # Elk_Day=Var_Mrand["Elk:Day",],
                                     # Day_Rep=Var_Mrand["Day:Rep",],
                                     # Elk_Rep=Var_Mrand["Elk:Rep",],
                                     # Rep=Var_Mrand["Rep",],
                   Residuals=Var_Mrand["Residuals",])


# All variance components (fixed and random) +++++
var_comp <- rbind(Var_Mfix, Var_Mrand)


# log of variance components +++++
log_var_comp <- t(log1p(var_comp)) # log(x+1)
log_var_comp <- cbind(id=1:dim(log_var_comp)[1],
                      log_var_comp)
log_var_comp <- as.data.frame(log_var_comp)
log_var_comp$id <- str_pad(1:dim(log_var_comp)[1], 2, pad = "0")
log_var_comp <- reshape2::melt(log_var_comp, id=c("id"))
log_var_comp$value <- as.numeric(log_var_comp$value)

names(log_var_comp) <- c("PC", "Effect", "Variance")
  

sum_var_comp <- rowSums(var_comp)

# Percent variance explained by each effect -----
var_comp_m1_abs <- sum_var_comp
var_comp_m1 <- var_comp_m1_abs*100/sum(var_comp_m1_abs)
names(var_comp_m1) <- names(var_comp_m1_abs)


# pdf(file.path(fig_path,"variance_components.pdf"),
    # height = 3, width = 3, pointsize = 12)

xlim <- c(0, ceiling(max(var_comp_m1)/10) * 10)
par(mar=c(0.5,2,3,0.5)) 
barplot(var_comp_m1, 
        main="Variance components \n percentage",
        xaxt="n",las=2, col=c(darkblue, turquoise,violetred ,
                              limegreen, gray67), border = NA,
        legend = names(var_comp_m1), 
        args.legend = list(x="topleft", inset=c(0,0),
                           box.lty=0,cex = 1,
                           y.intersp = 0.8))

# dev.off()

table_var <- rbind(var_comp_m1_abs, var_comp_m1)
rownames(table_var) <- c("Sum variance for all responses", 
                         "percentage of variation")
colnames(table_var) <- rownames(var_comp)
pander(table_var)
  Day Elk Day:Elk Residuals
Sum variance for all responses 75.09 175.7 113.9 438.2
percentage of variation 9.351 21.88 14.19 54.58
## Chapter graph

gg_color_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}

p <- ggplot(log_var_comp, aes(Effect, Variance, group = PC))+
  ggtitle("HSD - Log of variance components")+
  scale_color_manual(values = setNames(gg_color_hue(nPC),
                                       levels(log_var_comp$PC)))

p <- p + geom_point(aes(colour = PC))+
  geom_line(aes(colour = PC,linetype=PC),size=0.5)+ 
  ylim(c(0,5)) +
    ylab(label = "log(variance)") +
  theme(legend.text=element_text(size=10),
        legend.key.height=unit(0.7,"line"))+
   scale_linetype_manual(values = 1:nPC)+
  theme_classic()+
  theme(legend.key.width = unit(0.8,"cm"),
        plot.title = element_text(size = 17),
        axis.title.y = element_text(size = 17),
        axis.title.x = element_text(size = 17),
        axis.text = element_text(size = 13),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 17),
        legend.key.height = unit(0.5, "cm"))

# names(var_comp_m1)[5] <- "Residuals"
# p


# Thesis chapter output
# ggexport(p, filename = file.path(fig_path,
          # paste0(Vers,"HSD_variance_components.pdf")),
          # height = 5, width = 6.5, pointsize=20)
## Article journal graph
tab <- data.frame(Effect= names(var_comp_m1), 
                  pcvar = round(var_comp_m1,2))

var_comp_m1.table <- ggpubr::ggtexttable(tab,
                    cols = c("Effect", "Global var (%)"),
                    rows = NULL, theme = ggpubr::ttheme("classic",
                                                base_size = 10))

p <- p + annotation_custom(ggplotGrob(var_comp_m1.table),
                              xmin = 2.3, ymin = 3.8,
                              xmax = 0)
p

# ggexport(p, filename = file.path(fig_path,paste0(Vers,"ETS_variance_components_table.pdf")),
         # height = 6, width = 9, pointsize=16)

8.2 Method from Rousseau, 2011

fix_var <- (1/2)*2*(fixef_PC[2,]^2) # Method from Rousseau, 2011 for the fixed effects

var_comp <- rbind(fix_var, Var_Mrand)


sum_var_comp <- apply(var_comp, 1, sum)
pander( "sum of the variance components over the responses")

sum of the variance components over the responses

pander(sum_var_comp)
fix_var Residuals
111.4 438.2
sum_var_comp <- sum_var_comp*100/sum(sum_var_comp)
names(sum_var_comp) <- rownames(var_comp)
barplot(sum_var_comp)

pander(sum_var_comp)
fix_var Residuals
20.27 79.73

9 Bootstrap

# set up of the bootstrap
set.seed(2017)
nsim = 2000 # number of simulations

# name of the output 
name <- paste0("bootstrap_ELKtimeSeries_VSTT_repRandom.RData")

# Set up -------------------
# formulas without the effect to test
null_formulas <- list( 
  Volunteer = "~ Day + (1|Rep) ",
  Sampling =  "~ Day + (1|Elk)  ",
  # Tube = "~ Day + (1|Volunteer) + (1|Volunteer:Sampling)",
  Time = "~  Day + (1|Elk) + (1|Elk:Rep)")

null_effects <- names(null_formulas)

REML <- c(TRUE, TRUE, FALSE) # REML option for each null formula
names(REML) <- names(null_formulas)
# coding for the Volunteer:Sampling interaction
VS <- paste0("elk", design$Elk, ".", design$Rep) 
design$VS <- VS

9.1 Likelihood Ratio computation

##################################################
# True log-likelihood Ratio statistics      #
##################################################


# full model: MM_full
######################
# REML ++++++
loglik_PC_full_REML <- sapply(MM_full$merMod_obj, logLik, REML=T)
mean_loglik_PC_full_REML <- mean(loglik_PC_full_REML)

# ML (REML==FALSE) ++++++
loglik_PC_full_ML <- sapply(MM_full$merMod_obj, logLik, REML=F)
mean_loglik_PC_full_ML <- mean(loglik_PC_full_ML)



loglik_PC_full <- matrix(NA, ncol = nPC,
               nrow=length(null_formulas), byrow = TRUE)

for (i in 1:length(REML)){
  if (REML[i]==TRUE){
    loglik_PC_full[i,] <- loglik_PC_full_REML
  }else {loglik_PC_full[i,] <- loglik_PC_full_ML}
}


mean_loglik_PC_full <- rep(mean_loglik_PC_full_ML, 
                           length(null_formulas))

mean_loglik_PC_full[REML] <- mean_loglik_PC_full_REML



# Null models  
######################
res.parlmer_NULL <- vector("list", length = length(null_formulas))
names(res.parlmer_NULL) <- names(REML) <- names(null_formulas)


for (i in 1:length(null_formulas)) {
  # run parlmer
  res.parlmer_NULL[[i]] <- parlmer(design, 
                outcomes_step2, null_formulas[[i]], REML=REML[i])
}

# Save the results ------------
MM_PC_null <- lapply(res.parlmer_NULL, function(x) x[["merMod_obj"]])

ranNames <- sapply(res.parlmer_NULL, function(x) x[["ranNames"]])
Fixvarnames <- lapply(res.parlmer_NULL, function(x) x[["fixNames"]])

varcor_random <- vector(mode = "list", length = length(null_formulas))
fixef_PC <- vector(mode = "list", length = length(null_formulas))
modmat_fixed <- vector(mode = "list", length = length(null_formulas))
M0 <- vector(mode = "list", length = length(null_formulas))


for (i in 1:length(null_formulas)){
  varcor_random[[i]] <- sapply(MM_PC_null[[i]], 
                               function(x)
                                 as.data.frame(VarCorr(x))$vcov)
  
  rownames(varcor_random[[i]]) <- c(names(VarCorr(MM_PC_null[[i]][[1]])),
                                    "Residual")

  fixef_PC[[i]] <- sapply(MM_PC_null[[i]], fixef)

  modmat_fixed[[i]] <- model.matrix(MM_PC_null[[i]][[1]], type = "fixed")

  # intercept
  XM0 <- modmat_fixed[[i]]
  XM0[,-1] <- 0
  M0_PC <- XM0%*%fixef_PC[[i]]
  M0[[i]] <- M0_PC%*%t(spectra_PCA_loadings)
}


# Effect matrices computation ------------
index <- vector("list", length=length(null_formulas))
fixNames_int <- lapply(Fixvarnames, function(x) c("(Intercept)", x))

colnam <- lapply(modmat_fixed,  colnames)

index <- vector("list", length=length(Fixvarnames))
for (i in 1:length(null_formulas)){
  id <- c()
  for (k in 1:length(fixNames_int[[i]])){
    id <- c(id, grep(fixNames_int[[i]][[k]], colnam[[i]]))
  }
  index[[i]] <- id
}


Mfix <- Mfix_PC <- vector("list", length=length(null_formulas))
names(Mfix) <- names(fixNames_int)

for (i in 1:length(null_formulas)){
  XMfix = modmat_fixed[[i]]
  XMfix[,-c(index[[i]])] = 0
  Mfix_PC[[i]] = XMfix %*% fixef_PC[[i]] # Matrix of the Group effect
  # backtransform the PC to original coefficients
  Mfix[[i]] <- Mfix_PC[[i]]%*%t(spectra_PCA_loadings)
}

fitted_values_PC <- Mfix

######################
# compute the LRT 
######################

# objects initialisation ------------
Res_std_error_PC_null <- vector("list", length=length(null_formulas))
loglik_PC_null <- vector("list", length=length(null_formulas))
mean_loglik_PC_null <- c()
meanlog <- c()
sumlog <- c()




### compute the LRT ----------------------------
for (i in 1:length(null_formulas)){
  Res_std_error_PC_null[[i]]  <- sapply(MM_PC_null[[i]], sigma)

  # meanlog and sumlog
  loglik_PC_null[[i]] <- sapply(MM_PC_null[[i]], logLik, REML=REML[i])
  mean_loglik_PC_null[i] <- mean(loglik_PC_null[[i]])

  meanlog[i] <- -2*(mean_loglik_PC_null[i] /mean_loglik_PC_full[i])
  sumlog[i] <- 2*(sum(loglik_PC_full[i,] - loglik_PC_null[[i]]))
}

names(meanlog) <- names(null_formulas)
names(sumlog) <- names(null_formulas)


# Graphs ----------------------------

col1="blue"
col2="red"
par(mar=c(4,3,2,6))
dif <- vector(mode = "list")
for (i in 1:length(null_formulas)){
  # graphs
  mat <- cbind(loglik_PC_null[[i]], loglik_PC_full[i,])
  rownames(mat) <- paste0("PC", 1:nPC)
  col <- c(col1,col2)
  par(xpd=TRUE)
  barplot(t(mat), beside=T, ylab="Log-likelihood",
          cex.names=0.8, las=2, col=col, 
          main = paste("log-likelihoods",names(null_formulas)[i]),
          xpd=TRUE)

  legend("topright", legend = c("loglik restricted","loglik full"),
         fill = col, bty = "n", inset=c(-0.1,0))
  
   dif[[i]] <- 2*(mat[,2] - mat[,1]) # log-likelihood difference
}

difmat <- do.call(cbind, dif)
colnames(difmat) <- names(null_formulas)
col=c(darkblue,turquoise,violetred)
barplot(t(difmat), beside = TRUE, col=col,
        main = "Log-likelihood ratios")
legend("topright", legend = names(null_formulas), 
       title = "Removed effect", col = col, pch=15)

9.2 Booststrap tests

### Bootstrap ---------------------------------------

# bootstrapLT input arguments:
# MM_null = MM_PC_null$Sampling
# useREML=TRUE
# null_formula <- null_formulas[[2]]

 
bootstrapLT <- function(null_effect, useREML, MM_null, null_formula) {
  
  # simulate y from the null models --------
  nPC <- length(MM_null)
  simulatedY <- c()

  for (i in 1:nPC){
    ysim       <- unlist(simulate(MM_null[[i]], re.form=NA))
    simulatedY <- cbind(simulatedY, ysim)
  }
  

  y <- simulatedY
  dimnames(y) <- dimnames(outcomes_step2)

  # build restricted model -------
  # print("null")
  f_null <- parlmer(design, outcomes = y, null_formula, REML=useREML)

  
  # build full model -------
  # print("full")
  f_full <- parlmer(design, outcomes = y, form, REML=useREML)

  MM_f_null <- f_null$merMod_obj
  MM_f_full <- f_full$merMod_obj

  # LR -------
  loglikelihood_null <- sapply(MM_f_null, logLik, REML=useREML)
  mean_loglikelihood_null <- mean(loglikelihood_null)

  loglikelihood_full <- sapply(MM_f_full, logLik, REML=useREML)
  mean_loglikelihood_full <- mean(loglikelihood_full)

  sumlog <- 2*(sum(loglikelihood_full - loglikelihood_null))
  sumlog

}


sumlog_boot <- vector("list", length=length(null_formulas))

system.time(
for (i in 1:length(null_formulas)){

  null_effect <- null_effects[i]
  print(null_effect)
  sumlog_boot[[i]] = replicate(nsim,
                           bootstrapLT(null_effect = null_effect,
                           MM_null = MM_PC_null[[null_effect]],
                           useREML =REML[null_effect],
                           null_formula = null_formulas[[null_effect]]),
                           simplify = TRUE)

}
)

save(sumlog_boot, sumlog, file = file.path(data_path, name))
load(file=file.path(data_path, name))

9.3 Drawings

9.3.1 Histogram drawings

names(sumlog_boot) <- casefold(names(null_formulas), upper = FALSE)
# plot.new 
######################
# Time
######################
# pdf(file = file.path(fig_path, "ETS_hist_Time.pdf"),
    # width = 7, height = 5, pointsize = 15)

par(mar=c(2.1, 2.1, 3, 2), xpd=TRUE, mfrow=c(1,1))

df <- nPC *   1 # df for Time and Tube

m=hist(sumlog_boot$time, freq=F, breaks=100,
       xlab="Global Likelihood Ratio Statistic",
       xlim=range(sumlog["Time"], sumlog_boot$time),
       ylim = c(0,0.02),
       col = "gray75",border = "gray75",
       main = " Fixed Time effect", cex.main = 2.2)
lines(density(sumlog_boot$time), col=darkblue,lwd=2, lty=2)
lines(dchisq(seq(0,max(sumlog_boot$time)), df), col= "chartreuse4", lwd = 2)
points(sumlog["Time"], 0, col="red", pch=19, lwd=6)
legend("topright",
       legend = c(paste0("True GLRT: ", round(sumlog["Time"],2)),
                  "Kernel density", paste0("chi2 distrib. (df=", df,")")),
       col = c("red", darkblue, "chartreuse4"),
       lty=c(NA,2,1),pch=c(19,NA,NA),
       inset=c(0,0),box.lty=0, cex = 1.4, y.intersp = 0.8, lwd=c(4))

# dev.off()


######################
# Volunteer
######################
# pdf(file = file.path(fig_path, "ETS_hist_Elk.pdf"),
    # width = 7, height = 5, pointsize = 15)
par(mar=c(2.1, 2.1, 3, 2), xpd=TRUE, mfrow=c(1,1))
df <- nPC
m=hist(sumlog_boot$volunteer, freq=F, breaks=100,
       xlab="Global Likelihood Ratio Statistic",
       xlim=range(sumlog["Volunteer"], sumlog_boot$volunteer),
       col = "gray75",border = "gray75",
       main = " Fixed Elk effect",
       cex.main = 2.2)
lines(density(sumlog_boot$volunteer), col=darkblue,lwd=2, lty=2)
# lines(dchisq(seq(0,max(sumlog_boot$volunteer)), df), col= "chartreuse4", lwd = 2)
points(sumlog["Volunteer"], 0, col="red", pch=19, lwd=6)
legend("topright",
       legend = c(paste0("True GLRT: ", round(sumlog["Volunteer"],2)),
                  "Kernel density"),
       col = c("red", darkblue), lty=c(NA,2),pch=c(19,NA),
       inset=c(-0,0),box.lty=0, cex = 1.4, y.intersp = 0.8, lwd=c(4))

# dev.off()

######################
# Sampling
######################

# pdf(file = file.path(fig_path, "ETS_hist_Sampling.pdf"),
    # width = 7, height = 5, pointsize = 15)
par(mar=c(2.1, 2.1, 3, 2), xpd=TRUE, mfrow=c(1,1))
df <- nPC
m=hist(sumlog_boot$sampling, freq=F, breaks=100,
       xlab="Global Likelihood Ratio Statistic",
       xlim=range(sumlog["Sampling"], sumlog_boot$sampling),
       col = "gray75",border = "gray75",
       main = " Random Sampling effect",
       cex.main = 2.2)
lines(density(sumlog_boot$sampling), col=darkblue,lwd=2, lty=2)

points(sumlog["Sampling"], 0, col="red", pch=19, lwd=6)
legend("topright",
       legend = c(paste0("True GLRT: ", round(sumlog["Sampling"],2)),
                  "Kernel density"),
       col = c("red", darkblue), lty=c(NA,2),pch=c(19,NA),
       inset=c(0,0),box.lty=0, cex = 1.4, y.intersp = 0.8, lwd=c(4))

# dev.off()

9.3.2 bootstrap p-values

pval <- c()
for (i in 1:length(null_formulas)){
 pval[i] <- (sum(sumlog[i]<sumlog_boot[[i]])+1)/(nsim+1)
}

names(pval) <- names(null_formulas)
pander("p-values")

p-values

pander(pval)
Volunteer Sampling Time
0.0004998 0.1624 0.1959

9.4 p-value based on chi2 distribution for fixed effects

pander("true GLLR:")

true GLLR:

sumlog["Time"]
##     Time 
## 845.4103
df <- 1*nPC
curve(dchisq(x, df=df), col='red', main = "Chi-Square Density Graph",
          from=0,to=60)

pander("p-value Time")

p-value Time

pchisq(sumlog["Time"], df=df, lower.tail=FALSE)
##          Time 
## 1.782515e-151
pander("p-value Elk")

p-value Elk

pchisq(sumlog["Volunteer"], df=df, lower.tail=FALSE)
## Volunteer 
##         0

9.4.1 Graph

difmat <- do.call(cbind, dif)
difmat <- data.frame(PC=substr(rownames(difmat),3,4), difmat)
colnames(difmat) <- c("PC", names(null_formulas))
difmat <- gather(difmat, key=Effect, value = value, Time, Volunteer, Sampling)
difmat$Effect <- factor(difmat$Effect, levels = c("Time", "Volunteer", "Sampling"))
difmat$PC <- factor(difmat$PC, levels = as.character(c(1:40)))

class(difmat$PC)
## [1] "factor"
# difmat$Effect

LLRplot <- ggplot(data=difmat, aes(x=PC, y=value, fill = Effect)) +
  geom_bar(width=0.8,stat="identity", position=position_dodge(width=0.8)) +
  theme_classic()+labs(title="(Restricted) Log-likelihood Ratios") +  
  ylab(label="(R)LLR") +
  guides(fill=guide_legend(title="Removed effect"))
LLRplot

tab <- round(pval, 4)
tab[1:3] <- paste("<", round(pval, 4))
tab <- data.frame(Effect = names(tab), `p-value` = tab) #,
                  # chi2 = c("<5e-04","<5e-04","-", "-"))
pval_boot.table <- ggtexttable(tab,cols = c("Effect",
                                            "Boostrapped p-value"), 
                                            # "Chi2 test"),
                                 rows = NULL,
                               theme=ttheme(base_size = 9))

LLR_pval_plot <- LLRplot + annotation_custom(ggplotGrob(pval_boot.table),
                                                            xmin=33, xmax=33,
                                ymin=125, ymax=150) +
                                        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

# LLR_pval_plot <- ggarrange(LLRplot, pval_boot.table,
                # ncol = 1, nrow = 2,
                # heights = c(1,  0.3), common.legend=TRUE)

# ggexport(LLR_pval_plot, filename = file.path(fig_path,"ETS_LLR_pval_plot.pdf"),
         # height = 6, width = 5)

LLR_pval_plot

# journal article output

ggsave(LLR_pval_plot, filename = file.path(fig_path,paste0(Vers,"ETS_GLLR.pdf")),
         width=10, height=6, scale=0.85)
### Chapter graph

difmat <- do.call(cbind, dif)
difmat <- data.frame(PC=substr(rownames(difmat),3,4), difmat)
colnames(difmat) <- c("PC", names(null_formulas))
difmat <- gather(difmat, key=Effect, value = value, Volunteer, Sampling, Time)
difmat$Effect <- as.factor(difmat$Effect)
difmat$Effect <- factor(difmat$Effect,levels(difmat$Effect)[c(4,1,3,2)])

difmat2 <- cbind(Time=difmat[difmat$Effect=="Time", "value"],
                 # Tube=difmat[difmat$Effect=="Tube", "value"],
                 Volunteer=difmat[difmat$Effect=="Volunteer", "value"],
                 Sampling=difmat[difmat$Effect=="Sampling", "value"]
                 )

rownames(difmat2) <- c(1:dim(difmat2)[1])

gg_color_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}

# pdf(file.path(fig_path, paste0(Vers,"ETS_GLLR.pdf")),width=13,
    # height=8, pointsize = 20)



par(mar=c(4,4,4,1), cex=1.2)

# plotting settings -------------------------------------------------------
ylim <- range(mat)*c(1,1.5)
angle1 <- rep(c(45,45,135), length.out=4)
angle2 <- rep(c(45,135,135), length.out=4)
density1 <- seq(5,35,length.out=4)
density2 <- seq(5,35,length.out=4)
angle1 <- rev(angle1)
angle2 <- rev(angle2)
density1 <- rev(density1)
density2 <- rev(density2)

barplot(t(difmat2), beside=TRUE,col = gg_color_hue(4),ylab="(R)LLR", 
        xlab="PC", ylim=range(pretty(c(0,difmat2))),
         main = "(Restricted) Log-likelihood Ratios",
        angle=angle1, density=density1)
barplot(t(difmat2), beside=TRUE, add=TRUE, col = gg_color_hue(4),
        ylab="(R)LLR", xlab="PC", 
         main = "(Restricted) Log-likelihood Ratios",
        angle=angle2, density=density2)
legend("topright", legend = colnames(difmat2),
       title = "Removed effect:", col = gg_color_hue(4),
       fill=gg_color_hue(4),angle=angle1, density=density1,
       xpd=TRUE, inset=c(0,-0.1), bty="n")
par(bg="transparent")
legend("topright", legend = colnames(difmat2), 
       title = "Removed effect:", col = gg_color_hue(4),
       fill=gg_color_hue(4),angle=angle2, density=density2, bty="n",
       xpd=TRUE, inset=c(0,-0.1))

# dev.off()

# pdf(file.path(fig_path, paste0(Vers,"ETS_pval.pdf")))
pval_boot.table
# dev.off()

10 Session info

sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.3.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] phyloseq_1.38.0    RColorBrewer_1.1-2 car_3.0-12         carData_3.0-5     
##  [5] mdatools_0.12.0    tidyr_1.2.0        ggpubr_0.4.0       stringr_1.4.0     
##  [9] gridExtra_2.3      reshape2_1.4.4     cowplot_1.1.1      ggplot2_3.3.5     
## [13] pander_0.6.5       dplyr_1.0.8        plyr_1.8.6         lme4_1.1-28       
## [17] Matrix_1.4-0       MBXUCL_0.1         ezknitr_0.6       
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.4.1        igraph_1.2.11          splines_4.1.2         
##   [4] listenv_0.8.0          GenomeInfoDb_1.30.1    digest_0.6.29         
##   [7] foreach_1.5.2          htmltools_0.5.2        fansi_1.0.2           
##  [10] magrittr_2.0.2         cluster_2.1.2          recipes_0.1.17        
##  [13] globals_0.14.0         Biostrings_2.62.0      gower_1.0.0           
##  [16] bdsmatrix_1.3-4        prettyunits_1.1.1      colorspace_2.0-3      
##  [19] xfun_0.29              crayon_1.5.0           RCurl_1.98-1.6        
##  [22] jsonlite_1.8.0         survival_3.2-13        iterators_1.0.14      
##  [25] ape_5.6-1              glue_1.6.2             gtable_0.3.0          
##  [28] ipred_0.9-12           zlibbioc_1.40.0        XVector_0.34.0        
##  [31] genalg_0.2.0           spls_2.2-3             Rhdf5lib_1.16.0       
##  [34] future.apply_1.8.1     BiocGenerics_0.40.0    abind_1.4-5           
##  [37] scales_1.1.1           mvtnorm_1.1-3          DBI_1.1.2             
##  [40] plsVarSel_0.9.7        rstatix_0.7.0          Rcpp_1.0.8            
##  [43] progress_1.2.2         ropls_1.26.4           proxy_0.4-26          
##  [46] stats4_4.1.2           lava_1.6.10            prodlim_2019.11.13    
##  [49] htmlwidgets_1.5.4      ellipsis_0.3.2         farver_2.1.0          
##  [52] pkgconfig_2.0.3        nnet_7.3-16            sass_0.4.0            
##  [55] utf8_1.2.2             caret_6.0-90           labeling_0.4.2        
##  [58] tidyselect_1.1.2       rlang_1.0.1            munsell_0.5.0         
##  [61] phyclust_0.1-30        tools_4.1.2            cli_3.2.0             
##  [64] generics_0.1.2         ade4_1.7-18            pls_2.8-0             
##  [67] broom_0.7.12           evaluate_0.15          biomformat_1.22.0     
##  [70] fastmap_1.1.0          yaml_2.3.5             ModelMetrics_1.2.2.2  
##  [73] knitr_1.37             rgl_0.108.3            purrr_0.3.4           
##  [76] future_1.23.0          nlme_3.1-153           praznik_10.0.0        
##  [79] compiler_4.1.2         rstudioapi_0.13        MSQC_1.0.2            
##  [82] ggsignif_0.6.3         tibble_3.1.6           bslib_0.3.1           
##  [85] stringi_1.7.6          highr_0.9              lattice_0.20-45       
##  [88] nloptr_2.0.0           vegan_2.5-7            permute_0.9-7         
##  [91] multtest_2.50.0        vctrs_0.3.8            pillar_1.7.0          
##  [94] lifecycle_1.0.1        rhdf5filters_1.6.0     jquerylib_0.1.4       
##  [97] data.table_1.14.2      bitops_1.0-7           R6_2.5.1              
## [100] IRanges_2.28.0         parallelly_1.30.0      codetools_0.2-18      
## [103] boot_1.3-28            MASS_7.3-54            assertthat_0.2.1      
## [106] rhdf5_2.38.0           rprojroot_2.0.2        withr_2.4.3           
## [109] S4Vectors_0.32.3       GenomeInfoDbData_1.2.7 mgcv_1.8-38           
## [112] parallel_4.1.2         hms_1.1.1              rpart_4.1-15          
## [115] timeDate_3043.102      class_7.3-19           minqa_1.2.4           
## [118] clValid_0.7            rmarkdown_2.11         pROC_1.18.0           
## [121] Biobase_2.54.0         lubridate_1.8.0